
%***********************************************************************
% GMM.m
%***********************************************************************
% C.-Y. Cynthia Lin, 12 February 2010
%-----------------------------------------------------------------------
% The code provided is from the following paper:
%
% Lin, C.-Y. Cynthia.  (2013).  Strategic decision-making with information 
% and extraction externalities:  A structural model of the multi-stage 
% investment timing game in offshore petroleum production.  Review of 
% Economics and Statistics, 95 (5): 1601-1621.
%
% If the code is helpful to you in writing your own code, I would appreciate
% it if you would cite the above paper in your work and acknowledge me 
% for help with the code.  Thank you!
%-----------------------------------------------------------------------
%
% This function calculates the parameters of the oil production
% model using GMM by finding the set of parameters that best minimizes
% the moment conditions, given the desired estimator
% for Vce.  # moments = # params 
%
% =========================================================================
% ----------------
% Input Arguments:
% ----------------
%	lease_term
%		lease term
%	Te
% 		Last period of consideration for the exploration stage
% 		(this is the last period in which decisions can be made)
%	Td
% 		Last period of consideration for the development stage
% 		(all years from lease_term onwards are lumped together
% 		b/c tot_e_t, which depends on t, will not change from then
% 		onwards)
%	max_mkt_size
%		maximum market size (scalar)
%	beta
%		discount factor
%	tuple_mat
%	   	matrix of all the tuples of state variables that appear
%	   	in the data 
%	   (N_TUPLES X N_STATE_VARS)
%	n_per_tuple_vec
%	   	number of times ecah tuple of state variables was
%		reached in the data
%	   (N_TUPLES X 1)
%       tot_not_e_vec
%	   = total # who have not explored
%	     before year t when the state tuple at time t
%            is "tuple"	
%	   (N_TUPLES X Te_m)
%       tot_e_not_d_vec
%          = total # who are explored at t
%            (ie, by (t+1)) 
%   	     but have not developed
%	     before year t when the state tuple at time t
%	     is "tuple"	
%	   (N_TUPLES X Td_m)
%	n_d_vec
%	  = # who develop at time t 
%	    (ie, between t and (t+1)) when the state tuple 
%            at time t is "tuple" 
%	   (N_TUPLES X Td_m)
% 	Mn_emp
% 	 	empirical transition matrix from the point of 
% 		view of the owner of an unexplored tract 
%		that decides not to explore at 
% 		at time t.  The (i,j,t)th element gives the 
% 		empirical probability that the state
% 		tuple next period is tuple j, 
%		given that the tuple this period is
% 		tuple i and given that this tract owner 
%		decides not to explore this 
% 		period.
%	   (N_TUPLES X N_TUPLES X Te_m-1)	
%       Me_emp
% 		empirical transition matrix from the point of 
% 		view of the owner of an explored tract 
%		that decides not to develop at 
% 		at time t.  The (i,j,t)th element gives the 
% 		empirical probability that the state
% 		tuple next period is tuple j, 
%		given that the tuple this period is
% 		tuple i and given that this tract owner 
%		decides not to develop this 
% 		period.
%	   (N_TUPLES X N_TUPLES X Td_m)	
%       ge_emp
%		empirical probability of exploring 
% 		at time t given that the tract
% 		has not yet been explored and 
%		conditional on the publically
% 		available information at time t 
%		(but excluding the idiosyncratic
% 		private shock eps_i_t).
%	   (N_TUPLES X Te_m)	
%       gd_emp
%		empirical probability of developing
% 		at time t given that the tract
% 		has been explored but not yet developed 
%		and conditional on the publically
% 		available information at time t 
%		(but excluding the idiosyncratic
% 		private shock eps_i_t).
%	   (N_TUPLES X Td_m)	
%	avg_profitsifd_emp
%		average gross profits conditional on developing
%   		for each tuple reached in data
%	   (N_TUPLES X 1)
%       n_params
%		number of parameters
%       maxit_coarse    
%               number of iterations in optimization 
%               (coarse optimization)
%       tol_moms_coarse              
%		termination tolerance on moments in 
%		determining convergence
%               (coarse optimization)
%       tol_params_coarse              
%		termination tolerance on parameters in 
%		determining convergence
%               (coarse optimization)
%       tol_fp_coarse              
%		termination tolerance when determining the
%		fixed point for value functions 
%               (coarse optimization)
%	ini_diff_fp_coarse
% 		Initial difference between value functions before 
% 		solve for fixed point
%               (coarse optimization)
%       maxit_fine    
%               number of iterations in optimization
%               (fine optimization)
%       tol_moms_fine              
%		termination tolerance on moments in 
%		determining convergence
%               (fine optimization)
%       tol_params_fine              
%		termination tolerance on parameters in 
%		determining convergence
%               (fine optimization)
%       tol_fp_fine              
%		termination tolerance when determining the
%		fixed point for value functions 
%               (fine optimization)
%	ini_diff_fp_fine
% 		Initial difference between value functions before 
% 		solve for fixed point
%               (fine optimization)
%	Vce_estimator
%		one of: 'I', or 'G'
%	params0
%		vector of initial guesses for parameters 
%		(n_params X 1)
%
%
% -----------------
% Output Arguments:
% -----------------
%	sigma_mu
%		estimated parameter in exponential
%		distribution of mu
%	        (1 X 1)
%	sigma_eps
%		estimated parameter in exponential
%		distribution of epsilon
%	ce_factor
%		estimated coefficient of (ce_bin_t + 1) in exploration profit 
%		equation
%	        (1 X 1)
%	gamma_tote
%		estimated coefficient of tot_e_t in profit equation
%	        (1 X 1)
%	gamma_totd
%		estimated coefficient of tot_d_t in profit equation
%	        (1 X 1)
%	gamma_bid_bin
%		estimated coefficient of bid_bin_t in profit equation 
%	        (1 X 1)
%	gamma_cd_bin
%		estimated coefficient of cd_bin_t in profit equation 
%	        (1 X 1)
%	gamma_oilp_bin
%		estimated coefficient of oilp_bin_t in profit equation 
%	        (1 X 1)
%	gamma_constant
%		estimated constant term in profit equation 
%	        (1 X 1)
%	Vcn
%		estimated continuation value for each tuple
%               if decide not to explore 
%		an unexplored tract,
%		estimation using Vcn_estimator 
%	        (n_tuples X (Te_m))
%	Vce
%		estimated continuation value for each tuple
%		if decide not to develop
%		an explored tract,
% 		estimation using Vce_estimator 
%	        (n_tuples X (Td_m))
%	ge	
%		estimated exploration probability for each tuple
%	        (n_tuples X (Te_m))
%	gd	
%		estimated development probability for each tuple
%	        (n_tuples X (Td_m))
%       pi0_e
%		estimated deterministic component of profits
%               from exploring for each tuple,
%		estimation using desired estimator for Vce (and gd)
%	        (n_tuples X (Te_m))
%       pi0_d
%		estimated deterministic component of profits
%               from developing for each tuple
%	        (n_tuples X 1)
%	EVn
%		ex ante expected value at time 0 of an unexplored 
%		tract given the values of the tuples but not of mu,
%		estimation using desired estimator for Vcn (and ge)
%	        (n_tuples X 1)
%	EVe
%		ex ante expected value at time 0 of an explored but undeveloped
%		tract given the values of the tuples but not of epsilon,
%		estimation using desired estimator for Vce (and gd)
%	        (n_tuples X 1)
%	wtd_moments
%		weighted moments (1 X 1)
%
%***************************************************************************
function [sigma_mu, sigma_eps, ...
	  ce_factor, ...
	  gamma_tote, gamma_totd, ...
	  gamma_bid_bin, gamma_cd_bin, gamma_oilp_bin, gamma_constant, ...
	  Vcn, Vce, ...
	  ge, gd, ...
	  pi0_e, pi0_d, EVn, EVe, wtd_moments] = ... 
 	GMM(lease_term, Te, Td, max_mkt_size, ...
	  beta, tuple_mat, n_per_tuple_vec, ...
          tot_not_e_vec, tot_e_not_d_vec, ...
	  n_d_vec, ...
          Mn_emp, Me_emp, ge_emp, gd_emp, ...
          avg_profitsifd_emp, ...
	  n_params, ...
	  maxit_coarse, tol_moms_coarse, tol_params_coarse, ...
	  tol_fp_coarse, ini_diff_fp_coarse, ...
	  maxit_fine, tol_moms_fine, tol_params_fine, ...
	  tol_fp_fine, ini_diff_fp_fine, ...
	  Vce_estimator, params0);


% Global variables

    % Estimated continuation values
    global Vcn
    global Vce

    % Estimated investment probabilities
    global ge
    global gd
    
    % Estimated deterministic component of profits
    global pi0_e
    %global pi0_e_GG
    global pi0_d

    % Estimated ex ante values
    global EVn
    global EVe
    %global EVn_G
    %global EVe_G

    % Weighted moments
    global wtd_moments



% Determine the number of different tuples reached in the data
n_tuples = size(tuple_mat, 1);

% Determine the number of moment conditions
n_moments = n_params; 

% Initialize the weight matrix as an identity matrix 
% (n_moments x n_moments)
weight_matrix0 = eye(n_moments);

% Set the options for the coarse optimization
    options_coarse = optimset('Display', 'final', ...
                    'MaxIter', maxit_coarse, ...
                    'TolFun', tol_moms_coarse, ...
                    'TolX', tol_params_coarse, ...
                    'MaxFunEvals', maxit_coarse^2, ...
                    'LargeScale', 'off');

% Set the options for the fine optimization
    options_fine = optimset('Display', 'final', ...
                    'MaxIter', maxit_fine, ...
                    'TolFun', tol_moms_fine, ...
                    'TolX', tol_params_fine, ...
                    'MaxFunEvals', maxit_fine^2, ...
                    'LargeScale', 'off');

		    

% Find the values of the parameters that minimize the weighted moment
% conditions using the identity matrix as the initial weight matrix. 

    % Coarse optimization
    disp('--Coarse Opt--'); 
    params_ini_coarse = fminunc('wtd_moms', ... 
            params0, options_coarse, ...
	    	lease_term, Te, Td, ...
	        max_mkt_size, ...
		beta, tuple_mat, n_per_tuple_vec, ...
                tot_not_e_vec, tot_e_not_d_vec, ...
		n_d_vec, ...
          	Mn_emp, Me_emp, ge_emp, gd_emp, ...
                avg_profitsifd_emp, ...
   		weight_matrix0, tol_fp_coarse, ini_diff_fp_coarse, ...
		Vce_estimator);

    % Fine optimization
    disp('--Fine Opt--'); 
    params_ini_fine = fminunc('wtd_moms', ... 
            params_ini_coarse, options_fine, ...
	    	lease_term, Te, Td, ...
	        max_mkt_size, ...
		beta, tuple_mat, n_per_tuple_vec, ...
                tot_not_e_vec, tot_e_not_d_vec, ...
		n_d_vec, ...
          	Mn_emp, Me_emp, ge_emp, gd_emp, ...
                avg_profitsifd_emp, ...
   		weight_matrix0, tol_fp_fine, ini_diff_fp_fine, ...
		Vce_estimator);

% Let the params_ini_fine be the actual params, since we won't be using
% an optimal weight matrix
params = params_ini_fine;


% Separate params into its individual elements    
sigma_mu	= params(1,1);
sigma_eps	= params(2,1);
ce_factor	= params(3,1);
gamma_tote	= params(4,1);
gamma_totd	= params(5,1);
gamma_bid_bin   = params(6,1);
gamma_cd_bin    = params(7,1);
gamma_oilp_bin  = params(8,1);
gamma_constant  = params(9,1);


